home *** CD-ROM | disk | FTP | other *** search
/ PsL Monthly 1993 December / PSL Monthly Shareware CD-ROM (December 1993).iso / prgmming / dos / c / newmat.exe / NEWMAT5.CXX < prev    next >
C/C++ Source or Header  |  1991-07-30  |  8KB  |  285 lines
  1. //$$ newmat5.cxx         Transpose, evaluate etc
  2.  
  3. // Copyright (C) 1991: R B Davies and DSIR
  4.  
  5. #include "include.hxx"
  6.  
  7. #include "boolean.hxx"
  8.  
  9. typedef double real;                 // all floating point double
  10.  
  11. #include "newmat.hxx"
  12.  
  13. #include "newmatrc.hxx"
  14.  
  15. //#define REPORT { static ExeCounter ExeCount(__LINE__,5); ExeCount++; }
  16.  
  17. #define REPORT {}
  18.  
  19.  
  20. /************************ carry out operations ******************************/
  21.  
  22.  
  23. GeneralMatrix* GeneralMatrix::Transpose(MatrixType mt)
  24. {
  25.    if ((int)mt==0) mt = Type().t();           // type of transposed matrix
  26.    GeneralMatrix* gm1 = mt.New(ncols,nrows); gm1->ReleaseAndDelete();
  27.  
  28.    if (mt == Type().t())
  29.    {
  30.       REPORT
  31.       for (int i=0; i<ncols; i++)
  32.       {
  33.      MatrixRow mr(gm1, StoreOnExit+DirectPart, i);
  34.          MatrixCol mc(this, mr.Store(), LoadOnEntry, i);
  35.       }
  36.    }
  37.    else
  38.    {
  39.       REPORT
  40.       MatrixRow mr(this, LoadOnEntry);
  41.       MatrixCol mc(gm1, StoreOnExit+DirectPart);
  42.       int i = nrows;
  43.       while (i--) { mc.Copy(mr); mr.Next(); mc.Next(); }
  44.    }
  45.    tDelete(); return gm1;
  46. }
  47.  
  48. GeneralMatrix* SymmetricMatrix::Transpose(MatrixType mt)
  49. { REPORT  return Evaluate(mt); }
  50.  
  51.  
  52. GeneralMatrix* DiagonalMatrix::Transpose(MatrixType mt)
  53. { REPORT return Evaluate(mt); }
  54.  
  55. BOOL GeneralMatrix::IsZero() const
  56. {
  57.    REPORT
  58.    real* s=store; int i=storage;
  59.    while (i--) { if (*s++) return FALSE; }
  60.    return TRUE;
  61. }
  62.  
  63. GeneralMatrix* ColumnVector::Transpose(MatrixType mt)
  64. {
  65.    REPORT
  66.    GeneralMatrix* gmx = new RowVector; MatrixErrorNoSpace(gmx);
  67.    gmx->nrows = 1; gmx->ncols = gmx->storage = storage;
  68.    return BorrowStore(gmx,mt);
  69. }
  70.  
  71. GeneralMatrix* RowVector::Transpose(MatrixType mt)
  72. {
  73.    REPORT
  74.    GeneralMatrix* gmx = new ColumnVector; MatrixErrorNoSpace(gmx);
  75.    gmx->ncols = 1; gmx->nrows = gmx->storage = storage;
  76.    return BorrowStore(gmx,mt);
  77. }
  78.  
  79. GeneralMatrix* GeneralMatrix::Evaluate(MatrixType mt)
  80. {
  81.    if ((int)mt==0) { REPORT return this; }
  82.    if (mt==this->Type()) { REPORT return this; }
  83.    REPORT
  84.    GeneralMatrix* gmx = mt.New(nrows,ncols);
  85.    MatrixRow mr(this, LoadOnEntry); 
  86.    MatrixRow mrx(gmx, StoreOnExit+DirectPart);
  87.    int i=nrows;
  88.    while (i--) { mrx.Copy(mr); mrx.Next(); mr.Next(); }
  89.    tDelete(); gmx->ReleaseAndDelete(); return gmx;
  90. }  
  91.  
  92. GeneralMatrix* ConstMatrix::Evaluate(MatrixType mt)
  93. {
  94.    if ( ((int)mt==0 || mt==cgm->Type()) && cgm->Tag()==-1)
  95.       { REPORT return (GeneralMatrix*)cgm; }
  96.  
  97.    REPORT
  98.    GeneralMatrix* gmx = cgm->Type().New(cgm->Nrows(),cgm->Ncols());
  99.    MatrixRow mr((GeneralMatrix*)cgm, LoadOnEntry);//assume won't change this
  100.    MatrixRow mrx(gmx, StoreOnExit+DirectPart);
  101.    int i=cgm->Nrows();
  102.    while (i--) { mrx.Copy(mr); mrx.Next(); mr.Next(); }
  103.    gmx->ReleaseAndDelete(); return gmx; // no tDelete
  104. }  
  105.  
  106. GeneralMatrix* ShiftedMatrix::Evaluate(MatrixType mt)
  107. {
  108.    GeneralMatrix* gm=bm->Evaluate();
  109.    if ((int)mt==0)
  110.       { MatrixType mteqel(MatrixType::EqEl); mt = gm->Type()+mteqel; }
  111.    int nr=gm->Nrows(); int nc=gm->Ncols();
  112.    if (!(mt==gm->Type()))
  113.    {
  114.       REPORT
  115.       GeneralMatrix* gmx = mt.New(nr,nc);
  116.       MatrixRow mr(gm, LoadOnEntry); 
  117.       MatrixRow mrx(gmx, StoreOnExit+DirectPart);
  118.       while (nr--) { mrx.Add(mr,f); mrx.Next(); mr.Next(); }
  119.       gmx->ReleaseAndDelete(); gm->tDelete(); return gmx;
  120.    }
  121.    else if (gm->reuse()) { REPORT gm->Add(f); return gm; }
  122.    else
  123.    {
  124.       REPORT GeneralMatrix* gmy = gm->Type().New(nr,nc);
  125.       gmy->ReleaseAndDelete(); gmy->Add(gm,f); return gmy;
  126.    }
  127. }
  128.  
  129. GeneralMatrix* ScaledMatrix::Evaluate(MatrixType mt)
  130. {
  131.    GeneralMatrix* gm=bm->Evaluate();
  132.    int nr=gm->Nrows(); int nc=gm->Ncols();
  133.    if ((int)mt==0) mt = gm->Type();
  134.    if (mt==gm->Type())
  135.    {
  136.       if (gm->reuse()) { REPORT gm->Multiply(f); return gm; }
  137.       else
  138.       {
  139.          REPORT GeneralMatrix* gmx = gm->Type().New(nr,nc);
  140.          gmx->ReleaseAndDelete(); gmx->Multiply(gm,f); return gmx;
  141.       }
  142.    }
  143.    else
  144.    {
  145.       REPORT
  146.       GeneralMatrix* gmx = mt.New(nr,nc);
  147.       MatrixRow mr(gm, LoadOnEntry); 
  148.       MatrixRow mrx(gmx, StoreOnExit+DirectPart);
  149.       while (nr--) { mrx.Multiply(mr,f); mrx.Next(); mr.Next(); }
  150.       gmx->ReleaseAndDelete(); gm->tDelete(); return gmx;
  151.    }
  152. }
  153.  
  154. GeneralMatrix* NegatedMatrix::Evaluate(MatrixType mt)
  155. {
  156.    GeneralMatrix* gm=bm->Evaluate();
  157.    if ((int)mt==0) mt = gm->Type();
  158.    int nr=gm->Nrows(); int nc=gm->Ncols();
  159.    if (mt==gm->Type())
  160.    {
  161.       if (gm->reuse()) { REPORT gm->Negate(); return gm; }
  162.       else
  163.       {
  164.          REPORT
  165.          GeneralMatrix* gmx = gm->Type().New(nr,nc); gmx->ReleaseAndDelete();
  166.          gmx->Negate(gm); return gmx;
  167.       }
  168.    }
  169.    else
  170.    {
  171.       REPORT
  172.       GeneralMatrix* gmx = mt.New(nr,nc);
  173.       MatrixRow mr(gm, LoadOnEntry); 
  174.       MatrixRow mrx(gmx, StoreOnExit+DirectPart);
  175.       while (nr--) { mrx.Negate(mr); mrx.Next(); mr.Next(); }
  176.       gmx->ReleaseAndDelete(); gm->tDelete(); return gmx;
  177.    }
  178. }   
  179.  
  180. GeneralMatrix* TransposedMatrix::Evaluate(MatrixType mt)
  181. {
  182.    REPORT
  183.    GeneralMatrix* gm=bm->Evaluate();
  184.    if ((int)mt==0) mt = gm->Type().t();
  185.    GeneralMatrix* gmx=gm->Transpose(mt);
  186.    return gmx;
  187. }
  188.    
  189. GeneralMatrix* RowedMatrix::Evaluate(MatrixType mt)
  190. {
  191.    GeneralMatrix* gm = bm->Evaluate();
  192.    GeneralMatrix* gmx = new RowVector; MatrixErrorNoSpace(gmx);
  193.    gmx->nrows = 1; gmx->ncols = gmx->storage = gm->storage;
  194.    return gm->BorrowStore(gmx,mt);
  195. }
  196.  
  197. GeneralMatrix* ColedMatrix::Evaluate(MatrixType mt)
  198. {
  199.    GeneralMatrix* gm = bm->Evaluate();
  200.    GeneralMatrix* gmx = new ColumnVector; MatrixErrorNoSpace(gmx);
  201.    gmx->ncols = 1; gmx->nrows = gmx->storage = gm->storage;
  202.    return gm->BorrowStore(gmx,mt);
  203. }
  204.  
  205. GeneralMatrix* DiagedMatrix::Evaluate(MatrixType mt)
  206. {
  207.    GeneralMatrix* gm = bm->Evaluate();
  208.    GeneralMatrix* gmx = new DiagonalMatrix; MatrixErrorNoSpace(gmx);
  209.    gmx->nrows = gmx->ncols = gmx->storage = gm->storage;
  210.    return gm->BorrowStore(gmx,mt);
  211. }
  212.  
  213. GeneralMatrix* MatedMatrix::Evaluate(MatrixType mt)
  214. {
  215.    GeneralMatrix* gm = bm->Evaluate();
  216.    GeneralMatrix* gmx = new Matrix; MatrixErrorNoSpace(gmx);
  217.    gmx->nrows = nr; gmx->ncols = nc; gmx->storage = gm->storage;
  218.    if (nr*nc != gmx->storage)
  219.       MatrixError("Incompatible dimensions in CopyToMatrix");
  220.    return gm->BorrowStore(gmx,mt);
  221. }
  222.  
  223. GeneralMatrix* GetSubMatrix::Evaluate(MatrixType mt)
  224. {
  225.    REPORT
  226.    GeneralMatrix* gm = bm->Evaluate();
  227.    if (row_skip+row_number > gm->Nrows() || col_skip+col_number > gm->Ncols())
  228.       MatrixError("SubMatrix dimension error");
  229.    if (int(mt)==0) mt = MatrixType::Rect;
  230.    GeneralMatrix* gmx = mt.New(row_number, col_number);
  231.    int i = row_number;
  232.    MatrixRow mr(gm, LoadOnEntry, row_skip); 
  233.    MatrixRow mrx(gmx, StoreOnExit+DirectPart);
  234.    MatrixRowCol sub;
  235.    while (i--)
  236.    {
  237.       mr.SubRowCol(sub, col_skip, col_number);   // put values in sub
  238.       mrx.Copy(sub); mrx.Next(); mr.Next();
  239.    }
  240.    gmx->ReleaseAndDelete(); gm->tDelete(); return gmx;
  241. }   
  242.  
  243. void GeneralMatrix::Add(GeneralMatrix* gm1, real f)
  244. {
  245.    REPORT
  246.    register real* s1=gm1->store;
  247.    register real* s=store; register int i=storage;
  248.    while (i--) *s++ = *s1++ + f;
  249. }
  250.    
  251. void GeneralMatrix::Add(real f)
  252. {
  253.    REPORT
  254.    register real* s=store; register int i=storage; while (i--) *s++ += f;
  255. }
  256.    
  257. void GeneralMatrix::Negate(GeneralMatrix* gm1)
  258. {
  259.    // change sign of elements
  260.    REPORT
  261.    real* s1=gm1->store; real* s=store; register int i=storage;
  262.    while(i--) *s++ = -(*s1++);
  263. }
  264.    
  265. void GeneralMatrix::Negate()
  266. {
  267.    REPORT
  268.    real* s=store; register int i=storage; while(i--) { *s = -(*s); s++; } 
  269. }
  270.    
  271. void GeneralMatrix::Multiply(GeneralMatrix* gm1, real f)
  272. {
  273.    REPORT
  274.    register real* s1=gm1->store;
  275.    register real* s=store; register int i=storage;
  276.    while (i--) *s++ = *s1++ * f;
  277. }
  278.    
  279. void GeneralMatrix::Multiply(real f)
  280. {
  281.    REPORT
  282.    register real* s=store; register int i=storage; while (i--) *s++ *= f;
  283. }
  284.    
  285.